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Abstract 

Since the discovery of tumour initiating cells (TICs) in solid tumours, studies focussing on their role in 
cancer initiation and progression have abounded. The biological interrogation of these cells continues to 
yield volumes of information on their pro-tumourigenic behaviour, but actionable generalised conclusions 
have been scarce. Further, new information suggesting a dependence of tumour composition and growth 
on the microenvironment has yet to be studied theoretically. To address this point, we created a hybrid, 
discrete/continuous computational cellular automaton model of a generalised stem-cell driven tissue with 
a simple microenvironment. Using the model we explored the phenotypic traits inherent to the tumour 
initiating cells and the effect of the microenvironment on tissue growth. We identify the regions in 
phenotype parameter space where TICs are able to cause a disruption in homeostasis, leading to tissue 
overgrowth and tumour maintenance. As our parameters and model are non-specific, they could apply 
to any tissue TIC and do not assume specific genetic mutations. Targeting these phenotypic traits 
could represent a generalizable therapeutic strategy across cancer types. Further, we find that the 
microenvironmental variable does not strongly effect the outcomes, suggesting a need for direct feedback 
from the microenvironment onto stem-cell behaviour in future modelling endeavours. 

Author Summary 

In this paper, we present a mathematical/computational model of a tumour growing according to the 
canonical cancer stem-cell hypothesis with a simplified microenvironment. We explore the parameters 
of this model and find good agreement between our model and other theoretical models in terms of the 
intrinsic cellular parameters, which are difficult to study biologically. We find, however, disagreement be- 
tween novel biological data and our model in terms of the microenvironmental changes. We conclude that 
future theoretical models of stem-cell driven tumours must include specific feedback from the microen- 
vironment onto the individual cellular behavior. Further, we identify several cell intrinsic parameters 
which govern loss of homeostasis into a state of uncontrolled growth. 
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Introduction 

Heterogeneity among cancer cells within the same patient contributes to tumour growth and evolution. 
A subpopulation of tumour cells, called Tumour Initiating cells (TICs), or cancer stem cells, has recently 
been shown to be highly tumourigenic in xenograft models and have some properties of normal stem 
cells. Evidence continues to emerge that TICs can drive tumour growth and recurrence in many cancers, 
including, but not limited to, brain [1], breast [2] and colon [3]. These tumour types can be broadly classed 
as hierarchical as they have been posited to have a hierarchical organisation similar but not identical 
to non-neoplastic stem-cell (SC) driven tissues. In these hierarchical tumors, TICs can differentiate 
to produce non-TIC cancer cells or self-renew to promote tumor maintenance. As TICs have been 
demonstrated to be resistant to a wide variety of therapies including radiation and chemotherapy, the 
TIC hypothesis has important implications for patient treatments [4]. Specifically, the effect of current 
strategies on the tumor cell hierarchy should be defined, and TIC specific therapies are likely to provide 
strong benefit for cancer patients. 

In a simplified view of the tumour cell hierarchy, TICs can divide symmetrically or asymmetrically 
to produce two TIC daughters or a TIC daughter and a more differentiated progeny [5,6], respectively. 
More differentiated TIC progeny which still have the capability of cell division and are similar to transient 
amplifying cells (TACs) in the standard stem-cell model and are capable of several rounds of their own 
symmetric division before the amplified population then differentiates into terminally differentiated cells 
(TDs) which are incapable of further division. This mode of division and differentiation, which we will 
call the Hierarchical Model (HM) is schematized in Figure 1. 

In the HM, there are a number of cellular behaviours that govern the system. In this study, we choose 
to study three: the rate of symmetric versus asymmetric division of the stem cells (a), the number of 
rounds of amplification that transient amplifying cell can undergo before terminal differentiation (/?), 
and the relative lifespan of a terminally differentiated cell (7). While it is a simplification of reality to 
study only these three parameters and leave out others (for example: differing proliferation rates for 
the different cell types [7] or the differing metabolic demands of stem vs. non-stem cells [8]) rigorous 
quantification of these parameters has been extremely difficult to pin down experimentally and so the 
majority of the work to describe them has been in silico. Most germane to the loss of homeostasis is the 
work by Enderling et al. [9] which showed the changes to the size of a mutated tissue (tumour) as they 
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varied the number of rounds of amplification of TACs. Other recent work attempting to quantify the ratio 
of symmetric to asymmetric division in putative glioma stem cells was presented by Lathia et al. [10], who 
showed that this ratio can change depending on the presence or absence of growth factors, suggesting yet 
another method by which a tissue can lose or maintain homeostasis: in reaction to microenvironmental 
change. A critical limitation of in vivo lineage tracing performed to date is an inability to determine the 
impact of microenvironmental heterogeneity on TIC symmetric division. 

While the HM appears to be quite straight forward, there is growing evidence of complexity to be 
further incorporated into the model. There are likely to be differences in the extent of TIC maintenance 
or the ability of tumour cells to move toward a TIC state. TICs appear to reside in distinct niches 
suggesting there may be differences in the biology of these cells, but defining differences in TICs is 
limited by cell isolation and tumour initiation methods. Prospective isolation of TICs relies on surface 
markers, including CD133, CD151 and CD24 which can be transient in nature [11], due to modulation 
by the tumour microenvironment [12] or methods of isolation [13]. Characterisation of these sorted cells 
then requires functional assays including in vitro and in vivo limiting dilution assays [14]. 

As the importance of TICs becomes more evident as it pertains to aspects of tumour progression 
like heterogeneity [15], treatment resistance [16, 17], recurrence [18] and metastasis [19], the need for 
generalizable therapeutic strategies based on conserved motifs in these cells grows. We therefore aim 
to understand how the phenotypic traits discussed earlier (asymmetric division rate, allowed rounds of 
transient amplification and lifespan of terminally differentiated cells) and microenvironmental changes 
(modelled as differences in oxygen supply) effect resultant tissue growth characteristics. 

To this end, we present a minimal spatial, hybrid-discrete/continuous mathematical model of a hi- 
erarchical SC-driven tissue architecture which we have used to explore the intrinsic, phenotypic, factors 
involved in the growth of TIC-driven tumours. We consider parameters that involve the rates of division 
of the cells involved in the hierarchical cascade as well as micro-environmental factors including space 
and competition between cell types for oxygen. We present results suggesting that there are discrete 
regimes in the intrinsic cellular parameter space which allow for disparate growth characteristics of the 
resulting tumours, specifically: TICs which form tumours that are unsustainable, TICs that are capable 
of forming only small colonies (spheres), and TICs that are capable of forming fully invasive tumours in 
silico, just as we see diversity in biological experiments (Figure 2). 
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97 Results 

98 A systematic parameter exploration of the three key parameters relating to vascularisation of the domain, 

99 symmetric vs. asymmetric division (a) and progenitor division potential (/?) was performed. We also 

100 explored the parameter determining the lifespan of differentiated cells (7) and found that the only impact 

101 of longer lifespans is an increase in the amount of time before the simulations reach a steady state, but does 

102 not change the qualitative nature of the results. These results are summarised in Figure 3. Each of the 

103 three panels represents the results for a different degree of vascularisation (0.01, 0.05 and 0.1). A density 

104 of vascularisation of 0.05 would mean 12,500 oxygen sources in the domain. To determine the diffusion 

105 coefficient, we used the estimate of approximately 70 micrometers of effective oxygenation [20]. Each plot 

106 shows the total tissue size after 50,000 time steps as we change the proliferative potential of progenitor 

107 cells. Each of the lines shows a different ratio of symmetric vs asymmetric divisions. These results show 

108 that all these three parameters have a critical range where homeostasis is disrupted (tumourigenesis). 

109 Figure 4 shows examples of the typical results produced by this model. Although the proliferation 
no rates of all the cells remain the same, due to space constraints and the differences in a, the population 
in of TICs does not grow at the same rate as the non-stem population. Figure 4A shows an example of an 

112 unviable tissue (parameters: 0 = 0.001, a — 0.3, j3 = 50 and 7 = 1 day) where the vascularisation does 

113 not support the potential tissue size of that TIC, resulting in an area of hypoxia affecting the region that 

114 contains the TIC. That leads to the death of the stem cell and, eventually, the rest of the cells in the tissue. 

115 Figure 4B shows a case of slightly increased symmetric division, resulting in a dynamic homeostasis where 
lie cell birth and death is balanced so that tissue size remains relatively constant - which could represent 
117 the enigmatic dormant phase [9]. Finally, figure 4C shows an example where the system never achieves 
us true homeostasis. In this case a is slightly higher when compared with the previous example, suggesting 

119 a critical value at which overgrowth occurs. Over time, the number of TICs increases, allowing for the 

120 'tumour phenotype': unconstrained growth. Although this leads to areas of hypoxia, cells survive in the 

121 periphery of the blood vessels and keep growing until they take over the entire domain. A plot of cell 

122 number versus time for each of these three examples are plotted in Figure 5. 

123 Unsurprisingly, the higher the vascularisation of the domain the greater the tissue size it can support. 

124 Past a certain threshold, however, the difference becomes negligible and more remarkably, the qualitative 

125 dynamics are unchanged by any change in the microenvironment. The same effect is evident in the 
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126 other two parameters, the ratio of symmetric vs asymmetric division (a) of TICs and the proliferative 

12? potential of TACs (/?). Regardless of the vascularisation, disruption of homeostasis only occurs when 

128 the proliferative potential of TACs (/?) is below a maximum value of about 15. For values of symmetric 

129 division (a) above 0.3, the values for (3 in which this overgrowth occurs becomes even more restrictive 

130 with a range of approximately 10-15. 

131 Interestingly, we observed a conserved decrease in overall tissue size for the highest value of symmetric 

132 division, a = 0.5, when the progenitor cells were allowed only 5 divisions (/? = 5). We believe this 

133 phenomenon represents a situation where the tissue is not able to grow to its potential as the stem cells 

134 themselves occupy too much space, and never allow the progenitors to contribute as much as they could 

135 to the overall population. This is a supposition however, and deserves closer study. These results are 

136 summarized in Figure 3. 

13? Of note as well: in no simulation did we observe the 'tumour phenotype' for a value of a < 0.3, 

138 suggesting something akin to a 'phenotypic tumour suppressor' function for this parameter. As observed 

139 biologically [10], this ratio is highly susceptible to changes in microenvironment, suggesting an extension of 
uo this minimal model to include the microenvironmental factors measured in that study. How to incorporate 

141 the changes observed in that study into a mechanistic HCA model however, is not trivial, and we reserve 

142 it for a future extension of this work. Further, our current model exists only in two dimensions. 

143 While our quant iative parameters are based on experimentally derived values, the claims we 

144 make are largely qualitative abstractions, however, we stress that the specific quantitative 
us descriptions of cell fates are likely not yet accurate and could change if this model was in 
U6 three dimensions. 

147 Discussion 

us In this paper we have presented a simple two dimensional computational model of the HM of a TIC- 

149 driven tissue. Our results show that there are distinct regions in parameter space (that directly correlate 

150 to the intrinsic TIC phenotype space) that encode vastly different behaviour in the tissue (or tumour) 

151 arising from the TIC in question. These parameters represent different TIC phenotypes, and therefore 

152 do not represent any specific genetic mutation. In this way, we hope to generalise the intrinsic alterations 

153 which a TIC could undergo much in the same way that the hallmarks of cancer have generalised non 
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154 TIC-specific alterations [21]: our end goal being the identification of treatment strategies to target these 

155 phenotypes to slow or stop the progression of TIC-driven cancers. 

lse Because of the difficulties in understanding TIC specific traits in vivo, the biological data to support 

157 these conclusions remains sparse. There have been some carefully undertaken in vitro experiments on 

158 single TICs in glioblastoma, a highly invasive and malignant brain tumour, which suggest that TIC 

159 specific division behaviour (a in our model) is variable and changes based on environmental cues [10]. 

160 Further work has shown that the other microenvironmental cues, such as acidity [14] and hypoxia [22-28] 
lei can also alter the prevalence of the stem phenotype by utilising functional markers of sternness, but the 

162 mechanism for this increase is, as of yet, imperfectly understood. 

163 Our simulations do not show a significant TIC population dependence on vascular density (9), a 

164 surrogate for hypoxia, or a change in stem composition (see Supplementary Table 1), suggesting a flaw 

165 in the model. To rectify this, future iterations of this model should include direct feedback onto the 
lee cellular parameters from the microenvironment. We aim to parameterize this dependence by specific in 
16? vitro experiments designed to quantify this effect, rather than just elucidate its existence. Other future 

168 developments of this model should take into consideration the emerging body of work suggesting that the 

169 proportion of TICs within a tumour is directly affected by therapy and not just physiologic growth factor 

170 controls [29]. There is now evidence in several cancers to suggest that radiation increases the size of the 

171 TIC pool. Specifically, in breast cancer, it has been shown that radiation therapy induces non-stem cancer 

172 cells to de-differentiate into TICs [30]. Further, experimental studies have shown radiation increases the 

173 TIC pool in glioblastoma [31], which has often been attributed to radiation resistance associated with 

174 differences in cell survival [16]. A new study by Gao et al. [32], however, has shown in silico and in vitro 

175 that radiation can effect the symmetric to asymmetric division ratio (our intrinsic parameter a), yielding 

176 further clues about the mechanism of this TIC pool expansion. 

177 Dedifferentiation due to treatment related microenvironmental factors has not yet been considered in 

178 any spatial theoretical models. Dedifferentiation due to 'niche' specific factors was studied by Sottoriva 

179 et al. [33], who, using an agent based model, reported findings similar to ours: that the microenvironment 

180 made no significant change to the overall tumour growth dynamics. Beyond this single spatial study, the 

181 concept of SC dedifferentiation is gaining more and more attention in conceptual theoretical treatments 

182 [34] and has been modelled with a deterministic ordinary differential equation system for a well-mixed 

183 population of cells [35]. 
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184 We, as well as others, find that the HM of tissue growth does not completely capture all the necessary 

185 dynamics that characterise cancer growth - but there is still a great deal of understanding to be gained 

186 from studying this formalism. To this end, we have performed a study of the factors related to TICs 
is? driving this dynamic and have identified several key factors which promote increased growth of the 

188 resultant tumour. Motivated by Hanahan and Weinberg [21], who have simplified the myriad (epi)genetic 

189 alterations which a tumour can undergo into the hallmarks of cancer, we seek to distill the traits of 

190 TICs in a similar way. Specifically, our model suggests that the number of allowed divisions of TACs 

191 exhibits bounds outside of which tumour growth is unsustainable. This finding has been corroborated 

192 independently by recent work from Morton and colleagues [36]. Further, there is a specific balance of 

193 symmetric to asymmetric division which keeps tumours from overgrowing; almost acting as a phenotypic 

194 tumour suppressor. Indeed, changes in this ratio have been recently hypothesized to underlie the increasing 

195 stem pool in glioblastoma after irradiation [32], and could also hold a key to understanding tumour 

196 dormancy [9]. 

197 In summary, we have presented a minimal spatial Hybrid Cellular Automaton model of the HM of a 

198 TIC-driven tissue in which we have explored generalised TIC phenotypic traits and have identified several 

199 key cellular parameters which influence the overall tissue behaviour. While our model does capture a 

200 number of salient phenotypic characteristics of TICs that seem to be conserved, it fails to capture the 

201 recently observed changes in stem fraction secondary to microenvironmental perturbations. This is an 

202 indication that any computational model of a stem-hierarchical tissue, or tumour, built from this point 

203 on must not only include the physical microenvironment, but also feedback from the microenvironment 

204 onto the specific cellular parameters encoded in the HM. 

205 Therefore, this endeavour has identified the crucial point that the microenvironment must effect the 

206 behaviour of the cells within the HM, and also several conserved phenotypic hallmarks, which could be 

207 the result of any number of (epi) genetic alterations or microenvironmental perturbations. By focussing 

208 on mechanisms important for the HM of stem-cell driven tumour growth, we are seeking to identify 

209 common phenotypes which could be targeted in a variety of solid tumours in which TICs promote tumour 

210 maintenance - thereby reducing the number of therapeutic targets to a more tractable set. Only with 
2n this sort of distillation of the biological complexity inherent to cancer initiation (and indeed progression) 
212 can we hope to make progress against this disease. 
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213 Materials and Methods 

2u Our model is based on a hybrid, discrete-continuous cellular automaton model (HCA) of a hierarchically 

215 structured tissue. HCA models have been used to study cancer progression and evolutionary dynamics 

216 since they can integrate biological parameters and produce predictions affecting different spatial and time 

217 scales [15,33,37-40]. As shown in figure 6C, cells are modelled in a discrete fashion on a 500x500 2-D 

218 lattice. This comprises approximately 1cm 2 where we assume a cell diameter of 20 micrometers [41]. 

219 The domain has periodic boundary conditions but the simulations are stopped when a cell reaches one 

220 of the boundaries. Every time step, cells are iterated in a random fashion as to avoid any bias in the 

221 way that cells are chosen. Figure 6 A shows that, although all cells are assumed to have the same size 

222 and shape, they can only be one of three different phenotypes: TICs capable of infinite divisions, TACs 

223 which are capable of division into two daughters for a certain number (/3) of generations, and TDs which 

224 cannot divide but live and consume nutrients for a specified lifetime (7). Modes of division for TICs 

225 include asymmetric division (with probability 1 — a), which is division into one TIC daughter and one 

226 TAC daughter and symmetric division, which is division into two TIC daughters (probability a). 

227 The continuous portion of this model is made of up the distribution and consumption of nutrients (in 

228 this case modelled only as oxygen). Vessels, which are modelled as point sources and take up one lattice 

229 point (Vij in Equation 1), are placed randomly throughout the grid at the intiation of a given simulation, 

230 in a specified density (0). Each of these vessels supplies oxygen at a constant rate (A) which then diffuses 

231 into the surrounding tissue. The diffusion speed/distance is described by Equation 1, where 0(x,y,t) is 

232 the concentration of oxygen at a given time (t), and place (x, y), Dq is the diffusion coefficient of oxygen, 

233 A is the rate of oxygen production from a blood vessel, and \i t are the rates at which TIC, TAC 

234 and TD cells consume oxygen. The difference in time scales that govern the diffusion of nutrients and 

235 that at which cells operate is managed by updating the continuous part of the model 100 times per time 

236 step. During each update the oxygen tension in a given grid point is updated with the values of the 

237 surrounding cells using a von Neumann neighbourhood modulated by the diffusion rate (Do)- 



dO(x,y,t) 
dt 



= D 0 V 2 0(x, y, t) + \V x , y - jJLsS XiV - jJLpP x , y - HtT XiV 



(1) 



238 



Any simulation performed by this model can be characterised by the parameters found in Figure 7. 



239 



The most relevant parameters for the question we are trying to address are the following: 
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240 • Symmetric/asymmetric division rate of stem cells (a) 

241 • Vascular density of the tissue (9) 

242 • Number of allowed divisions of TACs (/3) 

243 • Lifespan of TDs (7) 

244 In each case, as can be seen in figure 6, a simulation is seeded with one TIC with a given set of intrinsic 

245 parameters (a, /?, 7) governing its and its progenys behaviour, which is placed in the centre of the 

246 computational domain. The domain is initialised with as many randomly placed oxygen source points 

247 (vasculature) as described by the vascular density parameter (0). 
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340 Figure Legends 




3 rounds of and death 

transient after Y 

amplification timesteps 



Figure 1. Cartoon representing the hierarchical model of stem-cell driven tissues. In this 
formulation, each stem can undergo two types of division, either symmetric (with probability a) or 
asymmetric (with probability 1 — a). Each subsequently generated transient amplifying cell (TAC) can 
then undergo a certain number (/3) of round of amplification before differentiating into a terminally 
differentiated cell (TD) which will live for a certain amount of time before dying (7 timesteps). It is 
these three parameters, which we assume are intrinsic to a given stem cell, which we explore in this 
paper. 
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Figure 2. Differential phenotypes in cultures enriched for brain tumour initiating cells. 

Bright field images of CD 133+ patient derived glioblastoma cell lines cultured in Neurobasal 
supplemented with EGF, FGF and B27, exhibiting striking phenotypic variability These differences 
highlight the heterogeneity present even in a highly controlled static environment between cells that are 
putatively the same. 

Tables 

Table 1. Supplementary Information. Raw stem and total cell numbers from several runs of the 
CA with varying parameter combinations. 
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Figure 3. Size of tissues vs. progenitor proliferative potential achieved by simulations 
using different levels of vascularisation and ratios of symmetric vs asymmetric divisions. 
Lines represent averages for each of the three realisations in each scenario. (Left). Low 
vascularisation density of 0.01 (Centre) Normal vascularisation density of 0.05 (Right) High 
vascularisation density of 0.1. In each of these cases, the maximum tissue size will depend on the right 
combination of a and /3. 



Downloaded from http://biorxiv.org/ on September 18, 2014 



15 




Figure 4. Three different examples of simulations resulting from the computational 
model. Each simulation represents one of the typical outcomes. Each begins with a single TIC 
seeded in the middle of the computational domain. In each situation the phenotype parameters are 
slightly different, resulting in (A) An unsustainable tissue (parameters: 6 = 0.001, a = 0.3, /3 = 50 and 
7 = 1 day), (B) A homeostatic tissue where the balance of stem cell sell renewal and progenitor 
proliferation leads to a tissue whose overall size remains relatively constant over time, possibly 
representing a dormant tumor (parameters: 6 = 0.05, a = 0.3, /3 = 15 and 7 = 1 day) and, (C) 
Neoplastic-like tissue where the tissue overgrows the computational domain (parameters: 6 = 0.05, 
a = 0.3, (3 = 5 and 7 = 1 day). On the right of the final time point, we have shown an 
example of the oxygen tension in the computational domain, n.b. these images are 
zoomed in to illustrate detail and do not represent the entire computational domain. 
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Figure 5. The three qualitiatively different tissue scale phenotypes plotted as cell numbers 
over time for the example simulations in figure 4. The black trace, representing the 
unsustainable simulation, grows quickly though never expands its stem population and then outstrips 
the available oxygen and collapses. The blue trace, representing the Homeostatic simulation, reaches a 
critical size and then maintains a steady birth-death balance. The red trace, representing the 
tumorigenic simulation, settles into an effectively linear trace on this log- log plot, suggesting power law 
growth. 
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Figure 6. Computational model description. (A) The model includes three different cell 
types: stem, progenitor and differentiated. All cell types interact with the micro environment in 
the form of oxygen tension. (B) The behaviour of each cell type is captured by a flowchart. The last 
segment with discontinuous arrows represents behaviour that is specific to the stem cells. (C) The cells 
are represented as agents inhabiting points in a grid in a 2D space with 500x500 grid points. Stem cells 
are represented as red points, progenitor as green and fully differentiated as blue. The vasculature is 
represented as oxygen source points in black. 



Parameter 


value 


Do (Oxygen diffusion) 


0.001728 


A (Rate of Oxygen production) 


1 


Ms, MP> M T 


0.0001 


a (Ratio of SC symmetric 
division) 


0.1, 0.3, 0.5 


P (TAC proliferative potential) 


1,5,10,11,12,13,14,15,16,1 
7,18,19,20,50,70,100 


y (Differentiated cell lifespan) 


1 


0 (Vascularisation) 


0.001, 0.01, 0.05, 0.1, 0.5 



Figure 7. Model parameters. 



